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Abstract 

We construct the phase diagram of a homogeneous two component Fermi gas with population 
imbalance under a Feshbach resonance. In particular, we study the physics and stability of the 
Larkin-Ovchinnikov phase. We show that this phase is stable over a much larger parameter range 
than what has been previously reported by other authors. 

PACS numbers: 03.75.Ss, 05.30.Fk, 34.90. +q 
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Trapped fermi gases offer us a wonderful opportunity to study strongly interacting fermion 
systems With Feshbach resonance, one is able to tune the effective interaction from 
weakly attractive for magnetic field above the resonance to strongly attractive below. This 
problem is also closely related to studies of quark and nuclear matter J2| 

Recent attention has shifted to systems with unequal populations E|. The latter 

problem is analogous to the physics of a superconductor under the influence of an external 
Zeeman field, which provides a chemical potential difference ^ — u\ = 2h between the two 
species | and j. In the weak-coupling limit, it was shown by Sarma [7j that the uniform state 
with population imbalance is unstable at zero temperature. By comparing the free energy of 
the normal state and the completely paired superconducting state, he concluded that, as the 
magnetic field is increased, there is a first order phase transition from the equal-population 
superconducting state (with gap A ) to the normal state at h = A /^/2. This implies that 
a system with given unequal numbers of the two spin species will either phase separate(with 
h = A /v^2), or be in the normal state (with h > Ao/\/2), depending on the imbalance. 
However, later, Fulde-Ferrell and Larkin-Ovchinnikov showed that there are better 
alternatives. Fulde and Ferrell (FF) considered a state with order parameter A(r) = Aoe* 9 ' r , 
that is, pairs with finite momentum. The order parameter has a spatially varying phase with 
however the magnitude still constant in space. They showed that in certain parameter space 
this state has lower free energy than the states considered by Sarma. Larkin and Ovchinnikov 
(LO), however demonstrated that, at least in the small order parameter limit, states with 
certain choices of sinusoidal variations of order parameter (such as A(r) = Aocos(gx) etc) are 
more energetically favorable than the FF state. Notice that in the LO state, the magnitude 
of the order parameter is no longer a constant in space. With decreasing magnetic field 
and hence decreasing population difference rid, many authors showed that the LO state 
evolves into a state with a set of domain walls. This has been demonstrated for one 10], 



two [ill ], as well as in three [12] dimensions and also for a d-wave superconductor [13|. In 
the small population imbalance limit, the order parameter has a constant magnitude almost 
everywhere in space (with value identical to the state with no population difference), except 
near the domain walls. The phase of the order parameter changes by 7r when these walls are 
crossed, and these are also the locations where the magnetization concentrates. This local 
magnetization arises from the occupation of bound states, available due to the presence of 
the domain walls. The physics of these domain walls is closely related to the 7r-j unctions in 
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SFS (S=superconductor, F=ferromagnet) junctions [yj], where for suitable parameters it is 
energetically favorable for the two superconductors to acquire a phase difference tt. 

Indeed, in three-dimensions, in the zero temperature and weak-coupling limit, Matsuo 
et al I2I demonstrated that the energy of the domain walls becomes negative for h above 
hd w = \ X 0.7457rT c ps 0.665A where A is the zero temperature gap for the completely- 
paired superfluid of equal populations. This hdw is less than the phase separation field 
h ps = Aq/v^ ~ 0.707A (where the free energy of the normal and the completely-paired 
superfluid states are equal). Hence the LO state must be more stable than both the uniform 
superfluid and the normal state (at least) for hdw < h < h ps . For h slightly above hdw, 
the domain walls are far apart. The system resembles the uniform paired state except for 
the occasional domain walls. This state thus allows for the possibility of arbitrarily small 
population difference between the two species. This problem is quite analogous to the lower 
critical field h c \ in type II superconductors, where the vortex energy becomes negative and 
vortices begin to penetrate the superconductor for h slightly above h c \. For increasing h, 
more and more domain walls are formed, the order parameter becomes more sinusoidal like 
and the state evolves smoothly to the picture given by LO 9J. 

For the resonant Fermi gas, it has been recognized that at intermediate coupling strengths, 
the uniform state with population imbalance must be unstable I2, 4]. Indeed this has also 
been demonstrated also by experiments 6|. However, the investigation into the actual phase 
diagram, that is, the question as to which phase appears where instead of the unstable 
uniform state, even for the case without a trap, cannot be regarded as complete, especially 
regarding the stability of the FFLO states. Many papers la, Il6j examined only phase 
separation, whereas some [19I considered in addition only the FF state with 

A oc e tq r . However, none of these works actually investigated the LO states. They conclude 
that the FF state exists only in a very narrow region next to the normal state in the weak- 
coupling regime, and deduce then that phase-separation occurs throughout the rest of the 
entire region where the uniform phase is unstable. In particular, they conclude that, for small 
population differences, phase separation occurs unless the dimensionless coupling constant 
1/kfd — ► —00. However, as seen already in the above paragraph, this is simply an artifact of 
the FF state, which cannot smoothly go into the completely paired equal population state 
[20! ] . We expect that at least, in the weak-coupling limit, the state should be LO for any 
population imbalance below that of the normal state. 



3 



In this work, we shall investigate the stability of the LO state for arbitrary strength of 
;he attractive interaction between the two species, thus generalizing previous works such as 
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12| beyond the BCS limit. We shall concentrate in particular on the small limit. 



More specifically, we compare the critical chemical potential difference h& w at which the 
domain wall energy becomes negative, to the critical field h ps for phase separation where 
the free energy of the normal phase becomes equal to that of the completely paired equal 
population superfluid phase. We find that for 1/kfa < (>) —0.845, hdw < i>)h ps . Therefore 
we conclude that, for small n<j, the LO state is more stable than the phase separated state for 
1/kfa < —0.845. By combining with previous results j^l, [5, [3] (and with some reasonable 
extrapolations), we sketch the appropriate phase diagram for our system. 
The mean-field Hamiltonian of our system can be written as 



H 



2m 



[A*(r)^V>T + cc] 



|A(f) 
9 



where a =|, J, for the two species, ip c {r) their corresponding field operators, A(r) a position 
dependent order parameter and g the coupling constant. For convenience below we shall 
also write /if = \i + h and /i; = // — h. 

This Hamiltonian can be diagonalized by the Bogoliubov transformation for a general 



inhomogeneous system 



2$: 




E 



Uj(r) v*j(r) 
-vj(f) u*j(r) 




(2) 



where aj, j3j are annihilation operators for quasiparticles with spin j and [ of the state 
labelled by a set of quantum numbers J. Uj(r), vj(r) obeys J d 3 f[\uj(r)\ 2 + \vj(r)\ 2 



1 



and the Bogoliubov-deGennes (B-dG) equation 




E, 




(3) 



Substituting eq (j2J) into eq §\§ and using eq (j3J) to eliminate the derivatives Vwj and 
Vf j's, we obtain the ground state free energy 

1 |A(f)| 2n 




-2EjvAr)' + 



V 2e„ 



m 



4:irh 2 a 



\A(r)\ 2 \+^(Ej-h)f(Ej~h) (4) 
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The last term arises from the occupation of states aj with Ej < h. Here / is the Fermi 
function, V the volume, and we have eliminated the interaction constant g in favor of the 
scattering length a via the relation i = ^J^) — y 277 wnere e J is the energy of the state 
J in the normal phase. (Strictly speaking we should label the quantum states in the normal 
phase by another set of quantum numbers J', but we shall not make such a distinction for 
simplicity in notations. See further below.) In eq (jlj), A(f) should be viewed as a variational 
parameter with respect to which F has to be minimized. 

We are interested in the domain wall energy for given /z, h and a. We evaluate this by 
calculating the energy difference between a state with a single planar domain wall and the 
uniform completely paired superfluid state. The latter is easy, since it is independent of h 
and the B-dG equation ([3]) can be solved by Fourier transform. We obtained, for wavevector 
k, the familiar quasiparticle energies E? = [(e^ — fi) 2 + lAj 2 ] 1 ^ 2 (where e& = 4^-) and hence 
the bulk free energy F s {fi) = J2d e k ~ H ~ E k + 7^p|A| 2 ] - i J^|A| 2 . Minimizing with 
respect to A gives the usual BCS gap equation — 47r ™ 2a A = A-^ 



1 in 

2E k h 2 k 2 



We can 



also compute the corresponding density via n = J 



d 3 k 



1 - 



and the corresponding 



Fermi wavevector kf = [Z^n) 1 ^ and express a in the dimensionless combination fc/a. These 
results are identical with those in 



22]. We had also used A(z) = const instead of a domain 
wall in the procedure described below and verified numerically that we indeed obtained the 
same free energy density for the uniform state. [Equating this free energy to that of the 
normal state at the same \x and h, we obtained the phase separation field h ps discussed.] 

Next we evaluate the free energy of a planar domain wall. For this, we put our system 
in a box of dimensions L x , L y = L x , L z = L. We assume that the order parameter varies 
only along the z direction and thus Fourier transform the x and y coordinates. Calling the 
resulting wavevector p, we rewrite (u j, vj) as 




u pJ {z) 



V p ,j{Z) 



(5) 



where f p is the component of r in the x-y plane, and j is now a quantum 
number for the z dependences. u p j(z),v Pt j(z) are dimensionless quantities obeying 
/ dz (\u pd (z)\ 2 + \v pJ (z)\ 2 ) /L = 1. Eq © becomes 

h 2 d 2 

(6) 
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where fl = // — The energies E p j and the wavefunctions u P: j(z), v p j(z) depend onp only 
through its magnitude p. The free energy can then be written as 

/I -i r I a / \ |2 "1 1 



(7) 

In this equation, we arrange the quasiparticle states for both the superfluid and the normal 
states as increasing function of the counting number j. We solve eq Q by discretizing the 
z coordinate. After the energies and the function v P j(z) are calculated, we put it in eq ([7]) 
and calculate the free energy. For A(z), we limit ourselves to the one-parameter ansatz 

A(z) = Aotanh f^j (8) 

with (3 as the variational parameter and = (2m/i) 1 / 2 . (Thus the width of the domain wall 
is given by f3k~ l .) Here Ao is the bulk order parameter for the given [i and a at h = 0, as the 
order parameter should approach this value far away from the domain wall. Our eq ([8]) is 



motivated by earlier investigations [ill, [12J, where their numerical results can be well fitted 
by a function of the form ([8]). (Since we are employing periodic boundary conditions in z, 
this ansatz actually introduces a sharp ((3 = 0) domain wall also at z = ±L/2. However, this 
can be taken care of easily by removing the contributions due to this extra domain wall.) 

We shall then input the ansatz eq ([8]) into eq ([6]) to solve for E p j. Our analysis of the free 
energy eq ([7]) is simplified by the following observations. At h — 0, since all quasiparticle 
energies are positive, f(E p j) = and the free energy is given simply by the integral in eq 
(J7j). The free energy at finite h of a given (3 is related to that of h — at the same /3 by 
simply adding the negative term ^2^j(E p j — h)f(E p j — h) due to the occupation of the 
quasiparticle states. 

An example for our results for = —1.0866, 1/kja = —1.0676 and A //i = 0.2 are 

discussed below. The bound state energies E b = E p j are illustrated in FigHJ We see that 



in general we have states below the continuum (E p j < |A | for /2 > 0, E Pi j < a/|/}| 2 + |A | 2 
for \x < 0. It is the bound states that are essential: as we shall see, the relevant values of h 
are below the gap edges). The (3 = results were checked against the analytical ones in the 
Appendix. The bound state energy for a given p decreases with the width of the domain 
walls, as one expects. The significance of this would be discussed again below. 
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The integrand of the first term in eq ((Tj) is plotted in Fig[2j (The double hump structure 
for small (3 0) is due to the fact that |A(z)| — > as z — ► 0: see eq (J7])). As expected, 
the integrand decreases to a constant corresponding to the free energy density in the bulk 
for distances sufficiently far away from the domain wall. We can then evaluate the domain 
wall energy Fd w at h = by simply integrating this excess contribution over z. At h = 0, 
Fdw is minimum at (3 ~ 3. Fdw is positive for all /3's, as expected since the uniform state 
should have a lower energy than a domain wall. For finite h, the free energies are evaluated 
by adding the negative term from bound state occupation as discussed above. The results 
are depicted in Fig [3j The free energy decreases due to the occupation of the bound states. 
Since the bound state energies are smaller for larger (3, the wall energy decreases faster for 
larger (3, shifting the f3 for minimum wall energy to larger values with increasing h (see 
FigH]). The domain wall energy becomes negative at sufficiently large h for all /3's. The 
important question is whether it will become negative for some (3 at a value of h which is 
less than that for phase separation. For the parameters under discussion, the domain wall 
energy first becomes negative for /3 ~ 5 at h = hdw ~ 0.1376//, which is less than the phase 
separation value h ps « 0.1391/i. We thus conclude that for 1/k/a = —1.0676, when rid is 
infinitesimally small, the system is in the LO state but not the phase separated state. 

At this point, it is of interest to compare our results more quantitatively with those in 12), 
who has assumed the quasiclassical limit in their calculations at the outset. They obtained 
the critical value for the sign change of the domain wall energy at h — hd w = \ X 0.7457rT c m 
0.665A (using A = 1.76T C ). Our hdw is close to their value if we simply rewrite their result 
as hdw/ '/•* = 0.665A /yU and substitute Aq//j, = 0.2. Also, their Fig 1 indicates a domain wall 
width near the critical field to be roughly given by 2A w//(7rT c ) 2 where Vf = kf/m is the 
Fermi velocity. Rewriting this expression as ~ 1.25 \jf~j (a^) ^/j 1 anc ^ simply subsituting 
the values of A //i etc we again find very good agreement. Thus, even though we did 
not assume the quasiclassical limit at the outset, for 1/kta = —1.067 our results can be 
understood from the quasiclassical limit calculations of 12| with extrapolations. 

For increasing 1/k/a, both hdw and h ps increase (and the (3 for the optimal domain wall 
decreases). However, hdw increases faster than h ps , and eventually hd w < h ps no longer 
holds. The situation is shown in Fig[5j If hd w > h ps , then the formation of domain wall is no 
longer favorable since phase separation already occurs at h slightly above h ps . We conclude 
then, for — > + , the transition from the LO state to the phase separation state occurs at 
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1/kfd ~ —0.845. This point is indicated by the point A in Fig [6j 

After locating this transition point for — > 0+, we now attempt to construct a phase 
diagram for general by combining the present results with previous ones in the literature. 
For n,i/n ^ and sufficiently negative 1/kfd, the system is in the normal state 

a a. 

We consider in turn the transition lines between this normal state and the LO and phase 



fl (c.f. 



separated states. Assuming that the transition to the LO state is second order 
below), the transition line between the normal and LO state can be found by solving the 
Cooper problem at finite wavevector q: 



"' 1 E 



Anh 2 a V 

k 



(9) 



^+e S -2/i h 2 k 2 

The transition line is determined by finding the optimal q corresponding to the weakest 
attractive interaction, i.e., the most negative 1/kfd. We note here that, since eq (Q is 
obtained by linearizing in the order parameter, the same equation determines the second 
order transition line into the FF state 9J. Our result can be turned into a line in our 
phase diagram (an equation of 1/kfd versus Ud/n) by using kf = (37r 2 n) 1//3 , with n = 
(2101 + h f/2 + { ^_ ^3/2] and m = j^^p This gives the long-dashed line 



in Fig Our numerical results (see also 231 ]) agree well with |5j. It should however be 



remarked that the transition from the normal state to the LO state has also claimed to 
be first order j^] at very low temperatures in three-dimensions. However, the difference 
between the actual transition line and the second order line is very small [241 ]. and we ignore 
this difference here. 

For the transition line between the normal state and phase separation, we equate [?| 
the free energy of the completely paired superfluid state Fs(n) to that of the normal state: 
F N (fx, h) = —V ^^i [(/i + /i) 5 / 2 + (/i — /i) 5//2 ] . This again yields a line 1/kja as a function 



of h/n and hence n^/n, shown as the solid line in figEJ Our numerical result (see also 23]) 
is also in good agreement with [5] . The two transition lines intersect at the point labelled by 
X at 1/kfd « —0.50. Therefore the transition is to the LO (phase separation) state if 1/kfd 
is less (larger) than this value. Interpolating between points A and X, we thus conclude 
that the system is in the LO state for the shaded region to the right of the line XA, whereas 
phase separation occurs for the shaded region to the left of this line. The transition lines 
between this phase-separated state and the homogeneous superfluid (Bose-Fermi mixture) 
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phase for positive 1/kja have already been discussed elsewhere in the literature ( {5], [15L Il6|. 
see also other references in and we would not repeat the details here. 

In conclusion, we studied the stability of the LO state in a homogeneous system, in 
particular in the limit of small population imbalance, by calculating the domain wall energies. 
Our investigation here is analogous to the determination of the lower critical field of the 
vortex phase of a type II superconductor. The determination of the detailed structures of 
the LO state, such as the question of the lattice structure of the domain walls (c.f. e.g. 24|, 
analogous to the vortex lattice structure in type II superconductors), as well as the phase 
diagram in the presence of a trap (studied already partially in 2a]). are left for the future. 

We thank C.-H. Cheng, C.-H. Pao and S.-T. Wu for their help in obtaining the transition 
lines from the normal state. This research was supported by the National Science Council 
of Taiwan under grant number NSC95-2112-M-054-MY3. 



Appendix — Bound states 

Due to the important role played by the bound states in the LO state, we give the analytic 
results for a sharp domain wall where A(z) = Aosgn(z) (corresponding to /3 — > 0). It should 
be noted that, even though |A(z)| is a constant throughout in space, bound states still exist 
due to the sudden change in phase factor of the order parameter from —71 to at z = 0. 
As already seen from eq (Q, the bound state energy is a function of /1 and p only through 
the combination ft. The B-dG equation ([6]) can be solved easily in this case since A(z) is 
piecewise constant. We require that the functions well as their derivatives 

be continuous at z = 0. The calculation is similar to solving the Schrodinger equation in 
piecewise constant potentials. It is straightforward and we simply state the results. It is 
convenient to divide them into two regions: 

(i) jx > 0. In this case states are bound when Ej, = E p j < |A |. We find that there is 
only one bound state (thus correspondingly one bound state for each p), with energy given 
by 

E t VTTW-i 



|A | 25 

where 5 = |A |//L We note here that the quasiclassical limit corresponds to |A | <C p, hence 
8 <S 1. In this limit, one can approximate the B-dG equation by the Andreev equation 26] 
and the bound state energy vanishes (irrespective of the detailed positional dependence of 
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A(z) so long as there is a it phase change: see e.g |27j and also 28]). For small | A n we 
have E b « |A | 2 / 2 /i- (cf. the mini-gap for bound states near the vortex core |2l|). With 
increasing |A | or equivalently decreasing p, ^/|A | increases. For 5 -> oo, E b -> |A |/\/2. 

(ii) /i < 0. In this case the continnuum starts at \/p 2 + | A | 2 . Again we find only one 
bound state, with 

E b VTTW+i 

|A Q | 25 1 j 

where 5 = |A |/|/2|. For /} = 0_ (5 -> oo), £ 6 -> |A |/V^ (cf. case (i) above). For 
decreasing /2 (decreasing S), E b increases. E b reaches |A | at 5 = 2. For p — > — oo, £J 6 
approaches the continuum from below. It is remarkable that bound state exists even for 
p < 0. We note that in this regime we have always E b > \p\. A similar situation occurs also 



for the bound states at a vortex 



29] for general coupling strength l/kja. 
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FIG. 1: (color online) Bound state energies for 1/kja = —1.0676, corresponding to = 
—1.0866, Ao//i = 0.2, for widths of the domain walls given by (3k~ l . Energies are normalized to n 
and the momentum p in the x-y plane normalized to ka- Inset: bound state energies in log plot. 
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FIG. 2: (color online) Integrand of eq j7|), in units of k^fiA where A is the area, as a function of 
position z (in units of k^ 1 ). Parameters are the same as in Fig [1] 
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FIG. 3: (color online) Domain wall free energy Fd w per unit area A (in units of k^fi) versus h (in 
units of Parameters are the same as in Fig[TJ Inset magnifies the region near F^ w = 0. 
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FIG. 4: (color online) Domain wall free energy Fd w per unit area A (in units of h&fi) versus (3. 
Parameters are the same as in Fig [TJ Inset shows the detailed behavior for h/fj, = 0.138. 
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FIG. 6: (color online) Phase diagram. N: normal state, BF: homogeneous superfluid (Bose-Fermi 
mixture) phase. Shaded regions to the right of XA: Larkin-Ovchinnikov state, shaded regions to 
the left of XA: phase separation. 
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